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FLAMES  IN  DUSTY  MIXTURES  -  THEIR 
STRUCTURE  AND  STABILITY 


J.  Buckmastet^ 

University  of  Illinois 
Urbana,  IL 

T.  Jackson^ 

Institute  for  Computer  Applications  in  Science  and  Engineering 
NASA  Langley  Research  Center 
Hampton,  VA 

ABSTRACT 

The  structure  and  stability  of  flames  in  dusty .  mixtures  is 
investigated.  The  presence  of  the  dust  leads  to  significant  transport  of 
energy  by  radiation  and  the  fimdamental  goal  of  the  analysis  is  to  explore  to 
what  extent  this  displaces  the  classical  non-hydrodynamical  stability 
boundaries  of  the  plane  deflagration.  An  approximate  description  of  the 
radiative  transport  permits  analysis  for  arbitrary  values  of  both  the  Planck 
length  and  the  Boltzman  number.  It  is  shown  that  the  pulsating  /traveling- 
wave  instability  usually  associated  with  values  of  Lewis  number  (Le)  bigger 
than  1  is  strongly  enhanced  by  the  presence  of  radiation  and  can  be  present 
even  if  Le  <  1.  On  the  other  hand  radiation  tends  to  suppress  the  cellular 
instability  normally  associated  with  values  of  Le  less  than  1.  The  latter  is 
consistent  with  preliminary  experimental  observations  of  Abbud-Madrid 
and  Ronney. 
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iDtroduciigp 


We  are  concerned  with  plane  flames  supported  by  mixtures 
that  contain  a  large  number  of  small  inert  particles  (dusty  gases). 
The  particles  substantially  influence  the  radiative  transport  and 
make  emission  and  absorption  energetically  significant.  Partial 
motivation  for  the  study  comes  from  recent  experimental 
observations  which  suggest  that  in  combustion  fields  of  this  kind 
cellular  instability  (the  Turing  instability  associated  with  values  of 
Lewis  number  less  than  1,  [1])  can  be  suppressed  [2].  A  discussion  of 
stability  and  the  effect  of  particle  loading  on  the  location  of  the 
stability  boundaries  in  the  wave-number  -  Lewis-number  plane  is 
the  central  contribution  of  this  paper. 

The  presence  of  particles  and  radiative  transport  presents 
serious  technical  obstacles  to  an  analytical  treatment,  and  these  must 
be  overcome  by  judicious  modeling.  We  start  by  assuming  that  the 
mass  loading  is  small  and  the  only  impact  of  the  particles  is  on  the 
radiative  transport.  The  bulk  particle  density  is(  4/3)  xd^  Np  (N  is 
the  number  density,  d  the  particle  diameter,  p  the  substantive 
particle  density)  and  the  radiation  effect  is  controlled  by  Nd^  so  that 
our  treatment  is  formally  valid  in  the  limit  N  eo,  d->  0,Nd^  fixed. 
(Then  0) 

A  strategy  for  dealing  with  radiative  transport  which  has  been 
usefully  employed  by  Joulin  and  his  colleagues  in  a  number  of 
important  studies  (e.g.  [3],-  [4])  is  to  assume  that  the  Boltzman 
number  is  small,  but  we  wish  to  avoid  the  restriction  here  since  it  is 
of  limited  applicability.  Instead  we  start  with  the  differential  - 
equation  approximation  (the  Eddington  approximation) 


a 


♦lvtV  (1) 

Since  the  nonlinearities  here  present  serious  difficulties  we 
approximate  the  Planck  length  L  by  a  constant  and  linearize  the 
emission  term,  replacing  Eqn.  (1)  by 

L^vfV-  q)  -  ^  +  18T?LVT  (2) 

where  Tc  is  a  constant  characteristic  (emitting)  temperature.  Then  q 
is  irrotational  so  that  we  may  introduce  a  generating  function  v  so 
that 

g-Vy.  (3) 

Equation  (2),  after  a  single  integration,  is  then  equivalent  to 

L^vSr  -  %  +  181?  L(T  -Tf)  (4) 

where  Tf  is  the  remote  cold-mixture  temperature. 

The  other  ingredients  of  our  model  are  familiar  and  of  proven 
merit.  They  are  embodied  in  the  equations 

p  (Yt  +  u Yx)  -  pD  -  Be  6(n) 

(5) 

pCp(Ti  +  uTx)  -XVT  -V V  +  QBe  5(n) 


Here  Y  is  the  mixture  mass  fraction,  T  the  temperature,  and  u  is  the 
steady  flame  speed  so  that  the  unperturbed  flame  is  fixed.  Reaction 
is  modeled  by  a  delta-function  whose  strength  is  an  Arrhenius 
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function  of  the  flame -temperature  T*;  n  is  the  distance  measured 
along  the  flame-sheet  normal  that  points  towards  the  burnt  gas.  And 
for  the  purposes  of  the  stability  analysis  we  adopt  the  constant- 
density  model  so  that  hydrodynamic  effects  are  discarded. 

Boundary  conditions  are: 

x->  T>Tf  .  Y-^Yf,v->0  (6) 

We  wish  to  construct  the  steady  planar  solution  of  the  system  (4)  - 
(6)  and  to  examine  the  linear  stability  of  this  solution  to  arbitrary 
perturbations.  To  this  end  it  is  convenient  to  introduce  a  flame-fixed 

coordinate  system  by  replacing  x  by  s  where 

s  s  X  -  Ry  . )  (7) 

and  the  flame-sheet  is  located  at  s  =  0.  Equations  (5)  for  s  0  are 
then 


p[Y,  +(u-F,)YJ  .pOvS"  , 
P  Cp[Tt  +  ( u  -F,)  TJ  - 


where 


S  (1  +FJ)  fds^+d^  /d/  -  a^y  /ayds  -Fyyd  /9S 


The  jump  conditions  at  the  flame-sheet  are: 


M  -  0  .[¥j  -  0 


4 


m  -  0  .[Y] .  0 


(10) 


(l  +Fy^[Tj  --(QBA)exp(-  BinT.)  ,(l  +Fy^  [YJ.(Bt>D)  exp(-  SRT.)  . 


The  steady  solution  is  constructed  by  seeking  solutions  for 


which 

F.O,  T4^(s).  Y.Y*(s)  .v->|^s) 
Consider  the  cubic 

aicoLf  -(<aLf  -(3a  +  P)(<dL)  +  3-0 

where  a  -X/LMCp  ,  P  «  IOT// MCpTc  . 


(11) 

(12) 

(13) 


a  is  the  familiar  length  X/M  Cp  divided  by  the  Planck  length;  p  is  a 
Boltzman  number. 

Equation  (12)  only  has  real  roots,  two  of  which  are  positive 
(ooi ,  a>2),  one  negative  ((03).  In  addition  we  define  (D4  b  y 


04  -pu /pD  s  M /pD  .  (14) 

Then,  without  using  the  connection  conditions  (10),  the  steady 
solution  has  the  form: 


^  O.  8  0,8 

s  <  0  T  —  Tf  +  Cl  0  +  C2  6 


—  Cl  |M  Cp  "Xcoijo  ■'C2(M  Cp  ■Xfl)2)® 


(15) 
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a>.  8 

Y»  -Yf  +  040  ^ 
s  >  0  T*  "Ta  +  Cse  ^ 

Vg  -  C3(m  Cp  -XQ>3)e“®“  (16) 

Y  =  0  . 

Here  Ta  s  Tf  +  CY  f /Cp  is  the  adiabatic  flame  temperature,  the 
temperature  at  s  -» «>«»  fixed  by  global  energy  conservation. 

The  various  constants  {Cj}  are  determined  by  use  of  the 

connection  conditions  (10  a-d)  whence 

C4  -  V,  (17) 

and,  with  the  definitions 

Zj-aLcoi  ,  Di  -Ci(Ta -Tf)  •  (18) 

we  have 

Di  +  D2  -  D3  ■  1 

Di  Zi  +  D2  Z2  -  D3  Z3  ■  1  (19) 

DiZi2+  D2Z22-D3Z32-  1 


Equations  (19)  have  solution 
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Di  -  (Z2  -  Wza  -  }iZ2  -Zi)  'Xza  -zi) 

□2  -  (Z1  -  )(Z3  -  )(Zi  -Z2)  ^23  -Zz) 

D3  -  -(Zl  -  )i(Z2  -  Xzi  -Z3)  •\Z2  -Z3)  (20) 

In  order  to  determine  the  mass  flux  M  (i.e.  the  flame  speed)  it  is 
convenient  to  deHne  Mq  by 

Mo  -  M(p  .()  (21) 

and  since  T.  >T^  at  this  limit,  Eqn.  (lOe)  is 

[TJ  -  (M„Cp/».)  (T, -TO  exp(E/2FT,  -  E/2R]  ,  (22) 

It  follows  that 

M/Mo  -  exp((Ei2RTa)D3[D3+Ta/(Ta-Tf)]i.  (23) 

D3  is  defined  by  the  roots  Zi  ,  Z2  and  Z3  which,  in  turn,  are  defined  by 
a  and  P,  both  of  which  depend  on  M.  Thus  Eqn.  (23)  is  an  implicit 
formula  for  the  flame-speed. 

The  limit  a  -*  (. 

It  is  of  interest  to  examine  the  realistic  limit  a  C,  when  the 
Planck  length  is  much  larger  than  the  diffusive  scale  (i.e.  the 
Reynolds  number  based  on  L  is  large).  Then  the  temperature 
distribution  is 


7 


SjsJI  T*.T(  +  (Ta  -Td[exp((MCp  /X)  ^  +  fjp  ,  4] , 

SJiSi  T*-T.  +  (T.-Td  f  .  (24) 

where 

f±(p , ^  +  ii  ‘''WII-I  p  ±  I  12) /  l]  sj 

On  the  scale  s  =  0(L)  there  is  a  radiative  preheat  and  postheat  zone 
separated  by  a  thin  region  (s  =  O  (X/MCp))  in  which  the  temperature 
increases  from 

T,  ^T,  +(Ta-Tf)p(p%li''  (25) 

to 

T.  ^Ta  +(Ta-T0p(p'+li'''^  (26) 

The  flame-speed  is  the  adiabatic  flame  speed  corresponding  to 
a  cold  mixture  of  temperature  Ti  so  that 

M/Mo  -  exp|(E(2RTa)p[p +Ta(Ta -Tf) .  (27) 

In  order  to  discuss  this  formula  it  is  convenient  to  introduce 
the  virtual  Boltzman  number 

Po  =  16TcVMoCpTc 


(28) 
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Using  CH4/air  adiabatic  flame  calculations  of  Giovangigli  and  Smooke 
[6]  and  with  Tc  equal  to  the  adiabatic  flame  temperature  we  can  plot 
variations  of  po  with  equivalence  ratio,  and  these  are  shown  in  Fig.  1. 
This  defines  a  realistic  range  of  values  of  Po  • 

Figure  2  shows  variations  of  M/Mo  and  P  with  Po  determined 
from  Eqns.  (27),  (28).  The  enhancement  in  the  burning  rate  was 
described  in  ref.  [3]  in  the  case  of  small  P  (p«0(RTa/8)  and  we  call  it 
the  Joulin  effect.  It  is  analogous  in  its  physical  origins  to  the 
enhancement  described  by  Weinberg  in  so-called  excess  enthalpy 
flames  [5].  Here  the  temperature  at  the  reaction  zone  is  raised  above 
the  adiabatic  flame  temperature  by  radiative  preheating,  whereas  in 
Weinberg’s  flames  the  increment  is  generated  by  providing  a  heat 
conduction  path  superior  to  that  of  the  gas. 

The  limit  a->0,j8->oo,aj3-O(1 

The  formulas  (24)  and  (27)  are  not  uniformly  valid  in  p. 
Indeed,  there  is  a  distinguished  limit  a  ^  0  ,p -» «»  ,aP  -  0(1)  In 
this  limit 


cdiL  ~  (2a)  '^[l  +Vl  +  4tp  ]  ,  q)2L  ~  &  /aP 

(29) 

(03L  ~  {2a)  [1  -Vi  +  4cp  ] 

Thus  the  radiative  preheat  zone  has  thickness  0(L^MCp  /^) 
whereas  the  postheat  zone  has  thickness  0(x/MCp).  The 
temperature  rises  slowly  from  Tf  to  Ta;  increases  rapidly  from  Ta  to 
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T.  -T,  +{T.-T()(1 +*p)  (30) 

and  then  decreases  rapidly  to  Ta.  The  burning  rate  is  given  by  the 
formula 

M  /  Mo  -  exp{(Ei2RTa)[l  +Ta(Ta  -Tf)  ’Vi  +4tp]  •’I  .  (31) 

It  is  well  known  that  in  the  limit  of  infinite  activation  energy 
the  reaction  zone  is  unstable  if  the  temperature  falls  off  towards  the 
burnt  gas  on  the  same  scale  (O  (X/MCp))  as  it  falls  off  towards  the 
fresh  gas  [7].  This  is  why,  for  example,  the  so-called  premixed -flame 
domain  of  diffusion  flame  burning,  for  which  the  burning  rate  is  a 
decreasing  function  of  Damkohle  number  (a  portion  of  the  middle 
branch  of  the  familiar  S-shaped  response),  is  unstable.  Thus  the 
steady  solution  described  in  this  section  may  not  be  physically 
realizable.* 

a  Not  Small 

When  a  is  not  small  we  return  to  Eqn.  (23)  in  order  to  calculate 
the  burning  rate.  Some  curves  are  shown  in  Fig.  2.  It  is  noteworthy 
that  for  high  particle  loadings,  when  a  is  not  small,  and  for  large 
values  of  the  virtual  Boltzman  number,  there  is  little  augmentation  of 
the  burning  velocity.  This  result  should  be  compared  with  the 

*  The  reaction  zone  instability  occurs  on  a  fast-time  scale  with  the  combustion 
Held  beyond  the  reaction  zone  unchanged.  It  is  not  captured  by  the  stability 
calculations  of  the  present  paper. 


saturation  observed  experimentally  [2].  Note  that  in  an  experiment 
Po  is  varied  by  changing  the  equivalence  ratio,  and  ao  can  be  varied 
with  Po  Hxed  by  changing  L  (i.e.  by  changing  the  particle  loading). 


Stability  Analysis 

We  seek  infinitesimal  perturbations  of  the  steady  solution  by 
writing 

T  -r  +eTo(s)e*‘^^"* 

Y  -  Y*  +eYo  (8)6"^^“^ 

V  -\|r®  +  eVo  (s)e  (32) 

_  *  iky  +0)1 

F  -  e$e 

where  e  ->  0. 

The  coefficient  functions  satisfy  the  equations: 
pCp  [to  To  +  u  To '  -8o)T*  ]  -X  [To  To  +  5k^  T®  ] 

-  [vo"''<^Yo  +  5k^V®] 

p[<dYo  +  uYo ’-8<dY«]  -pD[Yo"-k2Yo  +  8k2Y*]  (33) 

L2  [Y^"-k2y^  +  5k2y»]  -3Vj,+ IftTc^LTo 
Solutions  are: 

Th  8  Y-8  <D,  S  CD,8 

s  <  0  To  *Ai0  +A26  +fiGie  +5G20 
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Jt*  a 

Yo"C0  +8G70 


(34) 


Vo 

S2rQ  To 


n  *  «  ■'^2  *  ®2  * 

■■  Bi  0  +B20  *fSG4  0  '  +8650^ 

Y-8  Y^8  _  ®«8 

■A30  +A40  4-8630 


Yo=0.  (35) 

TTs  *  Ya  *  8 

Vg  ■iB3  0  ■fB40  -f8G6  0  ’ 

Hctc  >^2173  311^74  roots  of  the  quartic 

(yL)*  -a  'XlLp  -{iLf[2{kLf  +  3  -Mt  -koLu  -)]  +(YL)a  '\{kLf  +  ^ 

•►(kL)*+(kL)^(3 -fa  -f  (dLu  *))  +  3bt '^coLu  *''•0  (36) 

where  Re  (Y^,Y2)>t  and  Re  (73.7j<0  The  constants  Gi  through  G7 
are  known  (particular  solutions)  and  C,  Ai,  A2,  A3,  A4  are  unknown. 
The  constants  {Bi}  are  related  to  (Ai)  by 

B|  .  iaTc®L[L*(T|2-k2)  .  (37) 

Moreover  p  is  given  by  the  formula 

H  - 1  M  (pD)  '^  +  i  [m2(pD)  -2+  4{k2  +  pco  /po)]  (38) 

Thus  there  are  six  unknown  constants  (Ai  ,  A2  ,  A3  ,  A4  ,  C  ,8) . 

The  jump  conditions  (10),  when  linearized,  are: 
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[TJ  -  0  .[Vo]  -  0  .  [vo]  -  0  .  [Yo]  -  0  . 

(39) 

[To]  -  -  ^Ta  -Tf)To  /2RTf  aL  ,  [Yo]  -  B'f  MTo  /2RTf  pD  . 

These  provide  six  relations  amongst  the  six  constants  and  so  define  a 
linear  homogeneous  system  which  only  has  a  non-trivial  solution  for 
certain  values  of  the  eigenvalue  a>.  In  this  way,  stability  boundaries 
in  the  wave-number.  Lewis-number  plane  can  be  constructed. 

Representative  stability  results  are  shown  in  Fig.  3  for  an 
activation  energy  0  of  15.  Plotted  are  neutral  stability  boundaries 
for  (a)  a  =  1,  (b)  a  =  0.5,  (c)  a  =  0.1,  and  (d)  a  =  .02,  in  the  Lewis 
number  (Le«X/pDCp)  -  scaled  wave-number  (ak)  plane,  and  for 
various  values  of  the  Boltzman  number  (P). 

The  right  stability  boundary  (denoted  by  dashed  curves  in  Fig. 
3)  corresponds  to  a  pulsating  or  traveling  wave  instability  (Im  co  ^ 
0),  one  that  is  not  accessible  for  common  mixtures  when  p  =  0.  But 
radiative  effects  push  this  boundary  sharply  to  the  left,  to  commonly 
achieved  values  of  Le.  Figure  4a  shows  the  locus  of  the  intersection 
of  this  boundary  with  the  line  ak  =  c  in  the  Le-P  plane  for  c  = 
.001,  0.1,  0.5,  and  1.0,  when  a  =  1.0.  The  shift  to  the  left  (smaller  Le) 
with  increasing  P  is  eventually  reversed.  This  is  no  different  from 
the  effect  of  burner-induced  heat  losses  on  the  stability  boundary  in 
which  displacement  to  the  left  is  followed  in  due  course  by  reversal 
of  the  entire  curve  [8]. 

The  right  stability  boundary  is  not  shown  in  Fig.  3c  (a  =  .02). 
For  small  values  of  a  the  roots  are  extremely  sensitive  to  variations 
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in  P  and  there  are  serious  convergence  difficulties.  Sensitivity  to  a 
parameter  which  is  only  roughly  defined  (because  of  the  linearized 
transport  equation)  discouraged  us  from  examining  this  situation 
more  carefully.  If  it  is  to  be  done  it  should  be  done  in  the  context  of 
a  more  accurate  model  than  the  one  considered  here. 

The  analogy  with  burner-induced  losses  carries  over  to  the  left 
stability  boundary  (denoted  by  solid  curves  in  Fig.  3)  which 
corresponds  to  the  familiar  cellular  flames  and  for  which  o  -  0. 
There  is  a  strong  tendency  to  suppress  the  instability  by 
displacement  of  the  curve  away  from  Le  =  1.  Ronney  [2]  has 
reported  that  lean  methane/air  flames  are  smoothed  by  the  addition 
of  dust,  in  agreement  with  these  results.  Cellular  flames  in  dustless 
mixtures  are  usually  unsteady,  a  consequence  of  the  fact  that  the 
maximum  Lewis  number  on  the  stability  boundary  is  then  achieved 
when  k  s  0  (see,  for  example,  the  numerical  solutions  of  the 
Kuramoto-Sivashinsky  equation  in  [9]).  This  characteristic  does  not 
in  general  survive  the  addition  of  dust.  For  example,  when  a  =0.1,  P 
=  5  the  right  most  point  on  the  boundary  is  located  at  ak  0.42. 
Under  this  circumstance  we  anticipate  a  steady  cellular  structure  on 
a  scale  defined  by  this  special  value  of  the  wave  number.  The  same 
effect  is  seen  when  burner-induced  heat  losses  are  accounted  for  [8], 
as  seen  in  stationary  polyhedral  flames  [10],  [11]. 

Concluding  Remarks 

We  have  shown  that  the  introduction  of  significant  radiative 
transport  by  the  addition  of  dust  to  combustible  mixtures  can  have  a 
signiHcant  effect  on  the  thermal-diffusive  stability  boundaries.  The 
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pulsating  instability  can  be  strongly  enhanced  and  the  cellular 
instability  suppressed.  Suppression  of  cellular  instabilities  in  lean 
methane/air  mixtures  has  been  reported  by  Ronney  [2]. 

The  radiation  model  that  we  have  adopted  is  a  crude  one, 
although  it  appears  to  retain  the  essential  physics.  One  of  its 
characteristics  is  that  it  leads  naturally  to  the  definition  of  a 
Boltzman  number  that  contains  the  numerical  factor  16  rather  than  4 
(cf.  Eqn.  (13b)).  As  a  consequence  too  much  significance  should  not 
be  attributed  to  the  precise  values  of  p  that  appear  in  the  discussion 
and  in  the  figures.  The  order-of-magnitude  is  correct  however. 
Should  more  detailed  experimental  data  become  available  in  the 
future  it  might  be  appropriate  to  carry  out  a  more  accurate  analysis, 
together  with  more  detailed  exploration  of  the  stability  boundary 
locations. 
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Fig.  1  Virtual  Boltzman  number  vs.  Equivalence  Ratio  for  lean 
methane/air  flames,  calculated  using  the  flame  temperatures 
and  speeds  obtained  numerically  by  Giovangigli  and  Smooke  [6]. 


for  a  s  0, 0.1, 0.2, 0.5,  and  1.0  when  6  =  12,  TJT{  s  5.  Also  shown 
(dotted  line)  are  variations  of  the  Boltzman  ntimber  when  a  =  0. 


Locus  of  the  point  Re(ci))  -  0,  (xk  =  c  in  the  Lewis  number  -  p  plane. 
4a.  Right  branch,  0  =  15  ,  a  =  1.0,  Ta/Tf  =  5. 

4b.  Left  branch,  0  =  15  ,  a  =  1.0,  TaA'f  =  5. 


